% SWZmdd.M
% compute posterior mean, sd, and confidence interval (HPD)
% based on Metropolis-Hastings output

% Taeyoung Doh
% created: 02/02/2007
%--------------------------------------------------------------------1.0
% house cleaning
%--------------------------------------------------------------------
clear all;
close all;
diary off;

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec;
loaddata;
nobs = length(data); 

%--------------------------------------------------------------------
% prior specification
%--------------------------------------------------------------------
global prior_;
prior_ = feval(fnames_.prior,msel);

%lfile = [path_.para 'stat_mhrun' num2str(0) '.mat']; 
%load(lfile); 


global datasel psel 
% MHRUN
%  length of mhrun defines the number of columns in graph table
%	odd integer  (e.g., 1) for constant regime 
%	even integer (e.g., 2) for regime shifts 

mhrun    = 1; 
mhrunid  = num2str(mhrun','%2.0f');

block1   = 11;		% starting block
blockn   = 150;		% ending block
nsim     = 10000;	% number of simulation draws in each block

hpdprob  = 0.90;	% probability for highest posterior density

% for data density (modified harmonic mean)
hmax    = 20;

%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% parameter names
pnames = feval(fnames_.pnames,msel);


%--------------------------------------------------------------------
% means and standard deviations
%--------------------------------------------------------------------


% load path
 lpath = [path_.para 'mhrunpm' mhrunid '/' ];
%  lpath = [path_.para 'mhrun' mhrunid '/' ];  
% initialize output
ndraws     = 0;
drawmean   = 0;
drawsqmean = 0;
sumdrawsq  = 0;

t=1; 
draws     =[];
posts     =[]; 

% loop for block
for jblock=block1:blockn

	% load file
	lfile = [lpath 'mhdraw' num2str(jblock,'%04.0f')];
	load(lfile);

	% use only p.d. covariance matrix results for prior draws
	parasim = parasim(1:100:end,:);
    postsim = postsim(1:100:end); 
    
     
    draws = [draws;parasim];
    posts = [posts;postsim]; 
	% number of simulations in each block
	nsimul = size(parasim,1);

	% collect simulations
	drawblock = parasim;
    drawdim   = size(drawblock,2);
   % recmean(jblock-block1+1,:) = mean(totdraw);
	% compute sums of x(i) and x(i)^2 to be used to calculate means and s.d.
	drawmean   = drawmean + sum(drawblock,1);
	drawsqmean = drawsqmean + sum(drawblock.^2,1);
	sumdrawsq  = sumdrawsq + drawblock'*drawblock;
	ndraws     = ndraws + nsimul;
end

% Sims, Waggoner, Zha (2008) method 
[postmod,j]=max(posts);
paramod = draws(j,:);
covmod  = (draws-ones(size(draws,1),1)*paramod);
covmod  = cov(covmod); 
r       = zeros(size(draws,1),1); 
for i=1:size(draws,1)
r(i)  = (draws(i,:)-paramod)*inv(covmod)*(draws(i,:)-paramod)';
r(i)  = sqrt(r(i));
end 
inr   = find(r>0);
draws = draws(inr,:);
r     = r(inr); 
c1    = quantile(r,0.01);
c10   = quantile(r,0.1);
c90   = quantile(r,0.9);
v     = log(1/9)/(log(c10/c90));
a     = c1; 
b     = c90/(0.9^(1/v)); 
sigd  = cholcov(covmod); 
posts = posts(inr); 
crit  = quantile(posts,0.1); 
nmsim = 10000; % number of i.i.d draws simulated from the elliptical density to determine qL 
rr    = zeros(nmsim,1);
thetad= zeros(nmsim,length(paramod)); 
ind   = zeros(nmsim,1); 
qL =0;  

%while qL<0.0001
scale=0.25;   % 1 for other models, [0.1, 0.05, 0.15, 0.2, 0.25] for the four regime model with time varying trend 
sigd = scale*sigd; 
for k=1:nmsim 
u = rand;
rr(k) = (1-u)*a^v+u*b^v; 
rr(k) = rr(k)^(1/v); 
pl=-1000000;
while pl<=-1000000; 
thet    = mvnrnd(zeros(drawdim,1),eye(drawdim));
thetx   =  paramod+rr(k)*(sigd*thet')'/(sqrt(thet*thet')); 
[pp,pl]=rsfnmh(thetx',data);
end
thetad(k,:) = thetx; 
ind(k)  = (pp>crit);
end

qL=sum(ind)/length(ind);
disp('qL');
disp(qL); 
%end 
ns = size(draws,1);
 
hfac=120; 
mden =zeros(ns,1); 
for  k=1:ns
     
indd = posts(k)>crit; 
gden = log(gamma(0.5*drawdim))-0.5*drawdim*log(2*pi)-log(abs(det(sigd)));
gden = gden-(drawdim-1)*log(r(k))+log(v*r(k)^(v-1))-log(b^v-a^v);
hden = indd*exp(gden)/qL; 
mden(k) = (hden)/(exp(posts(k)-mean(posts)+hfac));
    
end
mmden = mean(mden); 
mdd2 = -log(mmden)+mean(posts)-hfac;
disp(        'Marginal Likelihodd : SWZ (2008)');
disp(  mdd2);  
save('SWZmddtrend.mat','mdd2'); 